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On  the  Computation  of  Simplicial  Approximations 
of  Implicitly  Defined  Two-dimensional  Manifolds  * 

Monica  L.  Brodzik  and  Wemer  C.  Rheinboldt 

Department  of  Mathematics  and  Statistics 
University  of  Pittsburgh,  Pittsburgh,  PA  15260 

Abstract:  A  method  is  presented  for  the  computation  of  a  simplicial  approxi¬ 
mation  covering  a  specified  subset  Afo  of  a  tivo-dimensinnaJ  manifold  M  in  R" 
defined  implicitly  as  the  solution  set  of  a  nonlinear  system  F(x)  —  0  of  n  —  2  equa¬ 
tions  in  n  unknowns.  The  given  subset  Mo  C  M  is  the  intersection  of  M  with 
some  polyhedral  domain  in  fj"  and  is  assumed  to  be  bounded  and  non-empty. 

The  method  represents  an  extension  (rf  a  local  simplicial  approximation  proce.ss 
developed  earlier  by  the  second  author. 

1.  Introduction 

For  nonlinear  mappings  F  :  R"  Jl"*,  n  =  m  -f  d,  d  >  2,  natural  conditions  exist 
that  guarantee  the  solution  set 

(1.1)  M  =  {i€/1";F(x)  =  0} 

to  have  the  structure  of  a  differentiable,  d-dimensional  manifold  M  in  fl".  We  present  here 
an  algorithm  for  computing  a  simplicial  approximation  that  covers  an  a  priori  specified 
subset  of  such  an  implicitly  defined,  two-dimensional  manifold. 

In  computer  aided  geometric  design  (CAGD)  and  related  applications,  manifolds  are 
often  parametrically  defined;  that  is.  as  the  image  set  {r  €  =  ♦(«,  u)}  of  some  kr-'wn 

parametrization  mapping  i  —>  R".  In  this  case,  various  triangulation  metliods  nave 
been  proposed,  (see  e.g.  [C93]  and  the  survey  (.B£192))  all  of  which  represent,  in  essence, 
extensions  of  techniques  developed  for  the  triangulation  of  flat  spaces. 

So  far  the  case  of  implicitly  defined  manifolds  (1.1)  has  not  received  as  much  attention. 
The  earliest  papers  appear  to  be  [A584j,  [A585],  and  [A<?87];  they  use  a  piecewise  linear, 
combinatorial  continuation  algorithm  to  construct  a  simplicial  complex  in  the  ambient 
space  fZ"  that  encloses  the  implicitly  given  d-dimensional  manifold.  The  barycenters  of 
appropriate  faces  of  the  enclosing  simplices  are  then  chosen  to  compute  a  global,  piecewise 
linear  approximation  to  the  maiufold.  However,  since,  in  general,  the  resulting  vertices  do 
not  lie  on  the  manifold,  this  does  not  represent  a  simplicial  approximation  in  the  standard 
seiue  of  combinatorial  topology. 

A  first  method  for  the  direct  computation  of  local  pieces  of  a  simplicial  approximation 
of  the  manifold  M  of  (1.1)  was  presented  in  (1187),  (fl88).  There  standardized  patches  of 
triangulations  of  the  tangent  spaces  T,  A/  of  M  are  projected  onto  the  manifold  by  smoothly 

*  This  work  was  supported  in  part  by  ONR-grant  N-00014-90-J-1025,  and  NSF-grant 
CCR-9203488 
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vHrying  projections  constr\icted  by  a  moving  frame  algorithm.  The  method  is  applicable 
to  manifolds  of  any  dimension  d  >  2  but  was  used  principally  in  the  case  d  =  2.  In  [/fo91| 
a  modified  version  of  this  method  for  the  case  d  2  was  considered.  It  uses  interpolating 
polynomials  to  predict  new  points;  but,  before  correcting  them  onto  the  manifold,  a  search 
is  conducted  to  determine  whether  the  new  points  may  cause  a  potential  overlap  with  any 
earlier  computed,  existing  triangle.  If  such  an  overlap  is  detected,  the  predicted  points  are 
identihed  with  appropriate  nearby  existing  vertices. 

Two  new  methods  were  developed  in  [AfAf93j  and  {//e93].  In  both  cases,  interest 
does  not  center  on  the  explicit  construction  of  a  simplicial  approximation  of  the  implicitly 
given,  2-dimensional  manifold  (1.1),  but  on  tesaellating  it  by  a  cell-complex.  In  [.ffe93] 
the  manifold  is  covered  by  overlapping  ellipsoidal  cells  each  of  which  is  obtained  as  the 
projection  of  a  suitable  ellipse  on  some  tangent  space.  In  [A/A/93J  a  complex  of  non- 
overlapping  cells  with  piecewise  curved  boundaries  is  constructed  by  tracing  a  fish-scale 
pattern  of  one<dimensional  paths  on  the  manifold.  Both  of  these  methods  appear  to  be 
intrinsically  designed  for  2-dimcnsional  manifolds. 

Here  we  present  an  extension  of  the  original  method  given  in  (HSS).  In  particular, 
the  process  is  globalized  to  allow  for  the  computation  of  a  simplicial  approximation  that 
covers  a  specified  domain  of  the  manifold.  The  algorithm  is  developed  for  the  case  d  =  2 
but  our  aim  was  to  use  tools  that,  in  principle,  can  be  generalized  to  higher  dimensional 
manifolds.  For  this  purpose,  the  mentioned  moving  frame  algorithm  is  replaced  by  a 
careful  consideration  of  the  orientation  of  the  triangles.  This  eliminates  the  cjjculation  of 
d-dimensional  singular  value  decompositions  which  can  become  costly  when  working  with 
manifolds  of  dimension  higher  than  two. 

In  the  original  algorithm  for  the  case  J  =  2,  the  patch  that  is  projected  from  the 
tangent  space  at  x  €  Af  onto  M  always  consists  of  a  hexagonal  neighborhood  of 
six  triangles  centered  at  x.  This  was  feasible  since  the  algorithm  was  only  applied  locally. 
But,  such  a  fixed  patch  is  likely  to  cause  local  overlaps  when  the  algorithm  is  applied  to 
larger  domains  of  M.  Thus  in  the  new  method  the  patches  are  constructed  adaptively  and 
are  allowed  to  have  fewer  than  six  triangles. 

The  algorithm  works  with  an  advancing  fremt  technique.  We  begin  with  a  point 
Xq  6  A/  and  add  it  to  the  database  that  stores  the  triangulation.  Then,  in  analogy  with 
the  original  method,  a  first  hexagonal  neighborhood  around  xq  is  constructed  in  the  tangent 
space  and  its  vertices  are  projected  onto  M  and  added  to  the  data  base.  The  starting 

front  of  the  process  is  formed  by  those  of  the  six  new  vertices  that  are  contained  in  the 
interior  of  the  given  domain  of  Af .  The  others  are  marked  as  exterior  vertices.  In  general, 
a  step  of  the  method  consists  in  the  selection  of  a  point  Xc  on  the  current  front  Then 
the  existing  triangles  incident  with  Xc  are  projected  onto  the  tangent  space  T«,A/.  This 
results  in  a  partial  neighborhood  of  x^  with  a  gap  that  still  has  to  be  closed.  If  the  gap 
is  too  small,  it  is  closed  by  identifying  its  two  open  edges,  otherwise,  it  is  divided  into  an 
optimal  number  of  triangles.  The  resulting  new  points  in  r,,Af  arc  then  projected  onto 
Af  and  added  to  the  data  base.  If  these  points  belong  to  the  given  domain  of  Af  they 
are  also  added  to  the  front,  otherwise  they  are  tagged  as  exterior  points.  Finally  the  step 
is  completed  by  removing  the  current  frontal  point  from  the  front,  since  it  now  has  a 
complete  neighborhood  of  triangles.  The  process  terminates  when  the  front  is  empty. 

The  resulting  simplicial  approximation  covers  the  given  domain.  It  is  not  difficult 
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to  adjust  all  exterior  points  onto  the  boundary  of  the  domain  although  this  may  resttlt 
in  needledike  triangles.  More  advantageous  is  to  use  a  local  Delaiiney-type  method  to 
effect  the  adjustment  onto  the  boundary  of  the  domain.  More  generally,  .such  a  Delaiiney 
approach  can  serve  also  to  improve  the  quality  of  the  entire  mesh.  We  shall  not  enter  into 
the  details  here.  As  noted  earlier,  our  atm  is  to  extend  the  method  to  implicitly  defined 
manifolds  of  dimension  larger  than  two.  This  is  the  topic  of  ongoing  work. 

In  Section  2  liclow  wo  give  a  brief  summary  of  some  barkgronml  material  needed  for 
the  development  of  the  method.  Then,  Section  3  outlines  the  data  structure  ust^l  here, 
and  Section  4  presents  a  detailed  description  of  the  algorithm.  Finally,  Section  5  shows 
several  numerical  examples. 

2.  Background 

Throughout  the  paper  F  :  V  C  R"  >-*  R’",  n  =  m  +  rf,  rf  =  2,  is  a  given  nonlinear 
mapping  of  class  C'  on  the  open,  connected  domain  T>  for  which  the  first  derivative  DF(j-) 
has  rank  m  for  all  x  6  D.  Then  is  is  well-known  that  the  solution  set 

(2.1)  M  =  {x€V;  F{t)  =  0], 

is  a  two-dimensional  C‘-manifold  in  R"  witlmut  boundary  (see,  e  g.,  |S79),  |/?S6]).  At  any 
I  €  A/  the  tangent  space  T,M  is  identified  with  kcrDF(x). 

Fbr  the  definition  of  the  subset  of  M  that  is  to  be  triangulated  wc  introduce  the 
hyperplanes 

(2.2)  =  {i€R";  Aj'(i-p»)  =  0},  ]t  =  l....ns, 

where  A,  €  is  a  imit  normal  vector  of  and  p,  €  R"  a  point  in  Hk.  Then,  with  the 
corresponding  half-spaces 

5t  =  {x6H";  Affi-ptl^O),  k=l,...n,, 

the  set  5  =  ns,  is  a  polyhedral  domain  in  R".  The  desired  subset  of  M  will  be  the 
intersection  Mo  =  SC\M.  We  assume  always  that  Mo  is  a  bounded  set  with  a  non-empty 
relative  interior.  Points  in  the  relative  interior  of  Mo  will  be  called  interior  points  of  Mo 
while  all  others  are  designated  as  exterior  points. 

As  in  [i{86|  we  introduce  at  any  “current”  point  x,,  g  A/  a  tangential  local  coordinate 
system.  For  this,  let  the  columns  of  I/'  g  /?“’'’  define  an  orthonormal  basis  of  the  tangent 
space  Tz,M  at  if  Then  the  implicit  function  theorem  applied  to  the  equation 

(2.3)  F(r, -tf/'y -b  DF(ic)^x)  =  0,  y  €  x  g  fi", 

guarantees  the  existence  of  open  neighborhoods  Ur  of  the  origin  of  R^  and  Vr  C  R"  of 
respectively,  such  that  for  any  y  €  ffc  there  exists  exactly  one  solution  z  of  (2.3)  with 
Xc  +  C^^y  +  DF{xe)^z  €  Vc  and  that  the  mapping  0  :  »-»  ^’(y)  =  Zy  is  of  cleiss  C* 

on  Uc-  Evidently,  wr  have  i/'(0)  =  0  and  /)0(O)  =  0  and 

♦  <t(y)  =  Zr  +  U^y  + DF(irfil>[y),  VygM,.  (2,4) 
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is  a  diffeomorphism  firom  We  onto  M  n  Vc.  In  other  words,  is  a  chart  of  M  at  jc  and 
we  call  4  a  tangential  local  coordinate  map  at  Xr- 

The  evaluation  of  ♦(p)  is  equivalent  to  projecting  the  part  x  =  (ic  +  U^y)  €  T,,Af 
orthogonal  to  T,^M  onto  M;  that  is,  to  solving  the  system 

F(*)  =  0 

(2.5)  i;‘^(i-(ic+y‘v))  =  o. 

Thus  in  general  we  require,  for  any  point  i  €  R"  sufficiently  close  to  Xc,  the  capability 
of  projecting  i  onto  M  orthogonally  to  T,^M.  There  are  several  ways  of  doing  this.  For 
example,  as  in  [R88|,  we  may  use  the  QR-factorization 

(2.6)  DF(x,f  =  (QuQi){^^ 

where  tgeQj  =  kerDF(Xc)^  and  set  =  Q-^.  Then  the  system  in  (2.5)  <*an  be  solved  by 
using  the  chord  Gauss-Newton  method 

x*+,  =  X,  -  DF{x,fiDF(x,)DF(x,f)-'F{xt).  k  =  1,2,... 

The  following  algorithm  incorporates  two  possibilities  for  providing  x,  namely,  (i)  x  6  R", 
X  ^  ic  is  directly  given,  or  (ii)  p  €  R*  is  provided  and  x  =  Xc  +  V‘y  is  still  to  be  computed. 
The  two  choices  are  passed  to  the  algorithm  via  the  input  variable  job,  with  job  =  1  for 
case  (i),  and  job  =  2  for  (ii). 

Proj:  Input:  Xj,  i,  p,  R,  Qi,  Qt,  job; 

if  ‘job  =  1’  then  i  ;=  x  else  i  ;=  Xt  +  Qjp  end 
while  'no  convergence’ 

solve  RT z  =  F(i)  for  z 

I  ;=  I  -  Qix 

Output:  X. 

For  all  X  sufficiently  near  ic  the  process  converges  to  a  unique  i*  6  M,  see  e.g.  [DH79). 

3.  Data  Structure  and  Problem  Definition 

Any  mesh  generation  algoritlun  depends  critically  on  the  data  structure,  for  which 
there  are  many  choices.  Since  we  restrict  ourselves  here  to  two-dimensioiml  manifolds, 
where  the  relationships  between  nodes  and  simplices  is  fairly  straightforward,  we  rho.se 
a  simple  data  structure  consisting  of  three  two-dimensional  arrays  xnod,  sim,  and  nod, 
for  the  node -coordinates,  the  simplex/node  incidences,  and  the  node/simplex  incidences, 
respectively.  Their  organization  is  summarized  in  Thble  3.1.  In  the  i-th  row  of  the  sim 
array,  the  indices  of  the  three  vertex  nodes  the  i-th  simplex  are  stored  in  an  order  that 
defines  the  simplex’s  orientation.  Any  two  coraecutive  simplices  listed  in  the  i-th  row  of 
the  nod  array  share  an  edge. 


Double 

Precision 

Integer 


Dimetmon 
n\axnod  rows 
nvnr  columns 


Contents  of  row  i 


Array 

xnod 


•im 


nod 


Integer 


mazsim  rows 
3  Columns 
maxnod  rows 
7  columns 


xnod(>,j)  —  j'*  coordinate  of  node  i, 

j  =  1, . .  ■  ,ntiar _ 

■iin(i,  j)  =s  index  of  j'*  vertex 

of  simplex  i,  j  =  1, . . . ,  3 _ 

nod(t,j)  =index  of  )'*  simplex 

incident  with  node  i,  j  =  .  ,k  <  6 

nod(t,  j)  =  0,  j  =  t  +  1, . . .  ,C 
nod(i,  7)  —  nodtyp  (See  later) _ 


Table  3.1;  Data  Structure 


The  operations  defined  on  the  database  are  as  follows: 

Addnod[f,  nodtyp]; 

Stores  coordinates  of  a  new  point  x  in  the  next  available  row  of  xnod,  and 
enters  the  point's  nodtyp  in  the  next  available  row  of  nod. 

Addsim[zi,  Z],  Zsj: 

Adds  a  new  simplex  to  the  aim  array,  and  updates  the  nod  array  by  adding 
the  simplex’s  index  to  the  rows  corresponding  to  the  vertices  zi.zj.xa. 

Equate[z,  Zi,  Zj,  zj; 

Identifies  two  computed  points  zi,zj,  which  are  incident  at  i,  by  replacing  zi 
with  the  projected  average  x  of  zi  and  zj,  removing  Zj,  and  then  updating  all 
three  tables  of  the  data  structure. 

Neighb[z]: 

Checks  if  the  given  point  z  is  connected  to  a  point  which  is  exterior  to  Mo- 

Note  that  in  this  data  structure  all  the  details  about  the  artiial  data  storage  and  ma¬ 
nipulation  had  to  be  included  in  the  software  package.  This  is  here  not  a  great  disadvan¬ 
tage  since  for  two-dimensional  manifolds  these  details  ate  fairly  simple,  and  our  resulting 
data  mauipultition  software  has  shown  to  be  acceptably  fast.  However,  when  generaliz 
ing  our  tiiangxilation  algorithm  to  higher-dimensional  manifolds,  we  will  use  the  relational 
database  management  system  SQL  to  keep  track  of  all  the  details  of  the  data  storage,  since 
the  relationships  between  nodes  and  simpUces  are  then  much  more  complicated. 

The  user  is  assumed  to  supply  the  following  three  subroutines  defining  the  problem: 

Fct[z,  F(z)|: 

Defines  the  function  F  in  (2.1)  and  returns  the  components  of 
F(x)  evaluated  at  the  given  point  z. 

Dfct(z,  DF(z)j: 

Defines  the  Jacobian  DF  of  F  and  returns  the  components  of 
DF(x)  evaluated  at  the  given  point  z. 
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Biida|x,  k,  bi,,  pt): 

Defines  the  hyperplanes  in  (2.2)  and  returns  the  components 
of  the  vectors  andpt  for  a  given  k  € 


Fbr  a  given  point  x  e  A",  the  following  algorithm  Cbkbnd  determines  whether  or 
not  X  belongs  to  the  polyhedral  domain  5.  In  addition,  it  computes  the  distance  dmin 
between  z  and  the  nearest  hyperplane  HkMtT- 


Chkbnd:  Input:  z,  number  of  hyperplanea  ns; 
for  k  =  l,ns 

{Pt>dt}  =  Bnds[jb]  /*Get  ps  €  Hs,  unit  vector  fcs  normal  to 
d  =  b^(z  -  pt)  /*Compute  signed  Euclidean  distance  d  from  z  to  fft.*/ 
if  d  <  0  then  /*z  does  belong  to  5:*/ 
knear  :=  k,  dmin  :=  d 

return 

else 

if  k  =  1  then 

knear  ;=  1,  dmin  :=  d  /’Initialize.*/ 
else  if  d  <  dmin  then 

knear  :=  k,  dmin  :=  d 

end 

end 

end 

Output:  knear,  dmin. 


4.  The  Triangulation  Algorithm 

The  triangulation  of  the  subset  Mo  =  5  n  M  of  the  manifold  M  begins  with  a  user- 
supplied  starting  point  zg  €  R"  which  need  not  be  on  M.  The  process  calculates  the 
QR-factorization  (2.6)  of  DF(xo)  and  uses  the  routine  Proj,  with  z  :=  zg  and  job  :=  1,  to 
project  Zg  unto  a  point  z„  €  Mo  ■  If  Proj  fails,  the  user  is  requested  to  supply  a  different 
starting  point.  Otherwise,  each  of  the  six  vertices 


(—  j).  (0.  H  (— .  j)' 


zIl\ 

V  2  ’  2 


(0.-H 


of  the  hexagonal  neighborhood  of  equilateral  triangles  around  the  origin  in  R^  is  mapped 
onto  the  affine  tangent  space  x„+T,„M  and  then  projected  onto  M,  using  again  Proj.  For 
h  either  a  user-suppUed  step  size  or  a  smaller  one  is  used  whichever  guarantees  successful 
projection  of  the  first  of  the  six  points  onto  M.  (This  first  successful  value  of  h  is  retained 
as  the  constant  step  size  throughout  the  remainder  of  the  triangulation).  The  projected 
points  inherit  the  connectivity  pattern  of  the  original  hexagonal  neighborhood  of  equilateral 
triangles  in 

Generally,  as  mentioned  earlier,  a  point  of  nod  is  either  an  interior  point  (of  Mg )  or 
an  exterior  point.  This  can  be  determined  by  means  of  Chkbnd.  Any  interior  point  is 
identified  as  a  frontal  point  if  it  does  not  yet  have  a  completed  simplicial  neighborhood. 
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Accordingly,  each  point  of  nod  is  classified  by  assigning  it  a  node  type  as  shown  in  Table 
4.1.  The  node  type  is  updated  as  the  triangulation  progresses. 


nodtyp 

Defiiu'Cio' 

-1 

FVontal  point  still  to  be  handled. 

U 

Interior  point  c<»mected  only  to  interior  points. 

1 

Interior  point  c<x)nected  to  an  exterior  point. 

2 

Exterior  point  connected  to  an  interior  point. 

Table  4.1:  Types  of  nodes  in  Af. 

The  frontal  points,  in  their  order  in  nod,  form  a  queue.  At  each  step  of  the  process 
the  topmost  point  is  removed  from  this  queue  and,  when  new  points  are  co.nputed  new 
frontal  points  may  be  added  at  the  end.  The  process  stops  when  the  queue  is  emi>ty. 

Let  Zc  denote  the  next  frontal  point  taken  from  the  queue.  Note  that,  after  the 
correction  of  the  starting  point,  r,  is  not  equal  to  Xo-  By  definition  a  portion  of  the 
simplicial  neighborhood  of  Xc  has  already  been  calculated  and  stored.  Figure  4.1  shows  a 
typical  example  of  this  existing  part,  and  Figure  4.2  gives  a  view  of  its  local  representation 
in  The  main  tasks  are  now  to  determine  the  gap  in  the  simplicial  neighborhood  that 
still  needs  to  be  closed,  to  decide  how  to  close  it,  and  finally  to  close  it. 


FVom  the  data  structure  the  two  simpliccs  ns*,,  and  n.Se„s  are  immediately  known. 
It  is  also  fiurly  easy  to  find  the  edge-defining  nodes  in  the  same  order  of 

consecutive  appearruice  around  as  the  incident  simplices  appear  in  nod.  Then  the  open 
edges  are  simply  ne*,,  :=  ji  and  ne,„j  ;=  ji+|.  However,  it  is  not  immediately  known 
whether  the  incident  simplices  n|,...,na  (with  :=  nod(i,m),  m  —  l,...,k)  appear 
clockwise,  as  in  Figure  4.1,  or  counterclockwise  around  x^  when  projected  back  onto  T,,  Af. 
Knowledge  of  this  ordering  around  Xc  is  essential  in  determining  whether  the  gap  is  the 
clockwise  angle  from  net,,  to  ne,„^  or  its  complement. 
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Central  to  the  determination  of  this  ordering;  is  the  unit  average  direction  vector  x 
defined  by 

1  **'  a 

The  gap  is  then  obtained  as  follofws.  For  x  and  the  normalized  open  edge  directions 


determine  the  local  coordinates  yi  ;=  U^Xi,  yj  :=  U^i,  and  yj  ;=  U^zt+i  in  f?*.  It 
is  expected  that  y^  will  point  in  a  direction  which  is  in  the  complement  of  the  gap  angle 
between  yi  and  y3.  (See  Figure  4.2).  The  reference  clockwise  rotation  angle  Orc/  of  yj 
into  Cl  :=  (1,0)^  is  calctilated  by  the  following  Givens-type  algorithm: 


Angle;  Input:  vector  (a,  k); 

nrm  :=  ||(a,  6)1|  /*  Euclidean  norm  of  (a, 6).*/ 

if  nrm  =  0  then  a  :=  0;  return  endif 
if  aba(a)  >  ais(6)  then 

if  6  =  0  then  ^  :=  0  else  ^  :=  arctan(b)  endif 
if  a  >  0  then 
a  :=  d 

if  k  <  0  then  a  :=  2x  —  ^  endif 

else 

a  ;=  X  —  d 

if  k  <  0  then  a  :=  x  +  endif 
endif 

else 

if  a  =  0  then 
^  :=0 


1 


9 


ebe 

^  ;=  arctan(f  ) 

if  a  <  0  then  ^  :=  endif 
endif 

If  4  >  0  then  a  ;=  f  -  0  ebe  a  :=  ^  +  ^  endif 

end 

Output:  a,  nrm. 

Let  A  denote  the  2x2  metrix  which  effects  a  clockwise  rotation  by  Or,/.  Then  A])i  and 
Ajfj  represent  clockwise  rotations  by  a,,/,  and  the  clockwise  rotation  angles  and  ^ 
of  the  resulting  vectors  dpi  and  dpi  into  ci,  respectively,  are,  effectively,  the  cl<M-kwise 
rotation  angles  from  pi  and  pi  into  pj.  (See  Figure  4.2).  If  ^i  <  ^i ,  then  the  simplices 
ni, . . .  ,nt  are  arranged  clockwise  from  n|  to  ns  around  the  origin  in  the  local  coordinate 
system.  Hence  the  gap  angle’s  magnitude  is  2r  —  /J| ,  and  its  orientation  indicator  l:or  is 
defined  to  be  1.  Otherwise  if  jd]  >  di>  then  the  simplices  arc  arranged  counterclockwise, 
the  gap  angle  is  0i,  and  we  set  kor  —  2.  The  following  algorithm  implements  this  gap 
determination  process: 

Agap:  Input:  center  point  z,  tangent  basis  U  ~  [ui,ui]  at  x; 

Order  z^, , . . , ,  x,,  /’Order  nodes  incident  at  x ,  so  x,, ,  z,^  define  open  edges.  •  / 
n«»«f  :=  111  ne,,s  ■=■  }p  /’Define  indices  of  open  edges.’/ 

Find  nss,,,  ns«nS  /’Get  indices  of  simplices  containing  open  edges.’/ 
i  =  2s=i(*>*  ~  *)»  *  ~  ®/ll*ll  /’Get  unit  average  direction  vector.*/ 

Is  :=  -  *)/llx),  -  x||,  k  =  l,p  /’Normalize  open  edge  directions.’/ 

y\'.=  U‘z\  /’Get  the  A’ local  coordinate  vector  of  i|.*/ 

Pi  :=  V^i  /’Get  the  A*  local  coordinate  vector  of  i.’/ 

Pj  ;=  f/^x,  /’Get  the  A*  local  coordinate  vector  of  Xp.*/ 
a  :=Angle[p]]  /’Get  clockwise  rotatiou  angle  a  of  pj  into  ei  :=  (1,0)^  .*/ 

^ _ /  cos  a  sin  a 

\  —sin  a  cos  Q 

ps  ;=  Aps,  i  =  1,2  /’Rotate  pi  and  pj  clockwise  by  a.’/ 
ds  :=Aiigle(psl  /’Get  clockwise  rotation  angle  0k  of  ps  into  ci,  t  =  1,2.*/ 
if  0]  <  01  then 

gap  :=  2r  —  0i  /’Define  gap  angle  which  needs  to  he  closed.*/ 
tor  :=  1  /’Define  orientation  kor  of  gap.  ’/ 

ebe 

gap  ;=  0t 
kor  :=  2 

end 

Output;  gap,  kor,  a,  ntk,f,  ne,„s.  nssep,  nSmW. 

Once  the  gap  angle  gap  and  its  corresponding  kor  value  have  been  determined,  the 
process  of  closing  the  gap  depends  on  the  magnitude  of  pap,  the  number  of  already  existing 
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MmpUces  rtti^  incident  at  ic,  and  a  fixed  minimum  acceptable  angle  size  amm.  If  the  gap 
is  small;  that  is,  if  gap  <  an,,,,  then  the  following  merge  algorithm  closes  the  gap  by 
“identifying”  the  open  edges  Ji  and  jj+i  and  then  replacing  the  nodes  defining  these  edges 
in  xnod  by  a  projection  of  their  average  ^{xj,  +  Xj,^, )  onto  M. 

Merge:  Input;  center  point  i;  other  endpoints  i|,  12  defining  open  edges; 

X  :=  l(n  +  *3) 

in  =iPro}[iyX,iiR,Qi,job  =  1]  /^Project  x  onto  A/.*/ 

Call  Equate|x,ri,r},7„i|  /^Replace  7i,X}  with  2,^  in  data  structure.*/ 

If  9^P  ^  Amin  and  riau  already  equals  the  maximum  allowable  value  of  6,  then  the  algo¬ 
rithm  stops  with  an  error  return.  Otherwise, 

n„,„  :=  mot  {  i  :  ^  >  a„i„,  1  <  «  <  (6  -  ) 

new  simplices  are  added  to  close  the  gap  at  i,  and  hence  to  complete  the  neighborhood 
of  simplices  around  Xe-  If  =  1,  then  only  one  simplex  is  added,  namely  the  one 

defiiirxl  by  the  already  existing  points  x,,  r,,,  and  The  indices  of  these  pmints  are 

entered  into  the  next  row  of  the  sim  array  in  an  order  such  that  the  orientation  of  the 
new  simplex  agrees  with  that  of  its  adjacent  simplex  If  n„e«.  >  1,  then  -  1 

new  points  are  needed  to  define  the  new  simplices  that  close  the  gap.  In  this  case,  the 
gap  -  angle  is  divided  into  nB,„  equal  angles,  and  the  new  points  in  the  local  coordinate 
system  are  defined  to  be  in  the  resulting  directions  and  at  a  distance  h  from  i,.  One  at  a 
time,  in  order  of  rotation  from  yj  into  the  gap,  these  points  in  Xc  +  T,,  Af  are  constructed, 
projected  onto  M  by  Proj,  and  added  to  the  database  by  Addncd.  As  each  new  point  on 
M  is  found,  Addaiin  adds  to  the  database  the  simplex  defined  by  the  new  point  and  the 
endpoints  of  the  open  edge  from  which  it  was  rotated.  At  the  last  new  point  x  €  Af  in  the 
gap,  a  second  simplex  is  added,  neunely  the  one  formed  by  i,  x„  and  x^, .  This  simplex 
completes  the  neighborhood  around  x^. 

For  some  new  direction  t,  Proj  may  fail  to  project  the  tangent  point  Xc  +  ht  onto  Af. 
In  this  case,  a  simple  continuation  process  is  started  along  the  direction  t  in  the  following 
way.  A  temporarily  smaller  stepsize  h  k  =  1,2,...  is  chosen  until  either  h  gets 

too  small  or  Proj  successfully  projects  the  corresponding  tangent  point  onto  Af.  If  A 
gets  to  be  smaller  than  some  minimum  ocreptable  stepsize,  the  algorithm  stops  with  an 
error  return.  Otherwise,  once  the  first  intermediate  point  xumr  G  Af  is  found,  further 
continuation  steps  are  taken  with  the  successful  stepsize  and  in  the  same  direction  t  until 
a  point  X  is  reached  where  the  sum  of  the  steps  exceeds  the  original  value  of  h.  This  x 
becomes  the  desired  new  point  to  be  added  to  the  database  together  with  the  simplex  it 
completes.  Then  h  is  reset  to  its  original  value. 

The  following  algorithm  summarizes  the  entire  triaiigulation  process. 

PITMAN:  Input:  Start  point  Xo,  suggested  step  A,  nunimum  gap  angle  amtn, 
total  number  of  bounding  hyperplanes  n*; 

DF(xof  =  (Oi.QzK?)  /‘Find  the  QR  decorap.  of  DF{xo).*/ 
x„  :=Proj[io,io,io,ff,Qi,Qz,l|  /’Project  x,  onto  Af.*/ 


n 


dmtn  :=  Chkbnd(zm.nt|  /'Check  if  x„  6  A/o>  i  e-  if  dniin  >  0.*/ 
if  'dtnin  <  0’  then  'error  return’  endif 
nod(l,  7)  :=  -1  /'Label  x„  aa  a  frontal  point.'/ 

Compute  the  neighborhood  of  six  simplices  of  x 
Remove  Xm  from  the  front, 
while  'FVont  is  non-empty' 

Get  the  next  frontal  point  x^. 

DF{x,)^  =  (<3i,Qi){J)  /'Find  the  QR  decomp,  of  DF{i,).V 
Call  Agap|ie,Q]|  /'Find  gap  and  open  edges  x,,  jj  at  ic  *l 
if  'gap  <  antin'  then 

Call  Merge[ic,xi,ii|  /'Identify  the  open  edges.'/ 

else 

if  ‘Xc  already  has  6  incident  simplices’  then 
‘error  return’. 

else 

Find  the  optimum  number  of  sectors  k  in  gap. 
if  ‘k  =  1’  then 

Call  Addsim(xr,X|,X'rj  /'Insert  one  simplex.'/ 

else 

Insert  k  simplices  into  the  gap. 

endif 

endif 

endif 

Remove  i,  from  the  list  of  frontal  points, 
eud  while 

Output:  Points  and  simplices  which  fonu  a  tciangulation  of  Mo. 


5.  Numerical  Experiments 

The  algorithm  PITMAN  described  in  the  last  section  has  been  implemented  in  Fortran 
77.  We  present  here  some  sample  problems  run  with  this  code. 

As  first  example  we  consider  the  intersection  of  the  unit  sphere  in  defined  implicitly 
by 

F(x)  =  x]  +  xl  +  xl  -  \  =  0. 


with  the  half  space 

5  =  {  X  €  R’  1  Is  >  -0.8  ). 

Figure  5.1  shows  a  view  of  the  triangulatiun  computed  by  PITMAN  using  the  stepsize 
h  =  0.3,  with  a  rotation  that  indicates  the  tnmeating  effect  of  the  hyperplaiie.  The 
algorithm’s  way  of  handling  local  overlap  causes  four  “seams”  of  elongated  triangles  on 
the  sphere.  This  could  be  remedied  by  a  Delaimay  improvement,  which  is  a  topic  of  ongoing 
work. 
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FIGURE  S.1 


The  second  example  arises  as  a  finite-eiement  model  of  the  deformation  of  a  thin, 
shallow  circular  arch.  This  test  problem  has  been  used  by  several  authors  and  can  be 
traced  back  to  (IVdQ].  We  use  the  same  model  formulation  as  in  (R86).  In  an  (r,tf)  polar 
coordinate  system  with  the  r-axis  as  vertical  direction,  the  unloaded  configuration  of  the 
arch  is  repicscntcd  by  the  segment  {  (r,S)  |  r  =  10,  — #o  <  #  <  #o  =  15“}.  Let  u  and  lu 
be  the  radial  and  axial  displacements.  For  pinned  ends,  the  dimensionless  total  potential 
energy  and  as.soriated  boundary  conditions  are  ipven  by 

j  ||(w'-u)  +  oo(«')’|  +oi(u'')’-ojpu 
u(tf)  =  U.(#)  =  u"{«)  =  0.  0  =  ±00, 

where  the  primes  denote  derivatives  with  respect  to  0.  Here  p  =  p{0)  is  the  dimensionless 
radial  load  and  oq,  O)  ,  O]  are  dimensionless  constants. 

As  in  [H8G]  wc  use  the  finite  element  approximation  consisting  of  a  uniform  mesh  with 
eight  elements.  The  problem  was  run  with  the  following  load  function 

{p  -  4p(«'  -  0)/0o,  V  V  ~  JOo  Si  9  <  v\ 

ft  +  4p(i'  -  9)/9(»,  if  v  <  0  <  V  +  j9o  ; 

0  otherwise. 


considered  already  in  [R88|,  where  v  and  ft  are  control  parameters.  In  other  words,  the 
load  is  a  piecewise  linear  hat  function  which  has  the  value  ft  ot  0  =  v  and  is  zero  outside 
the  interval  of  width  O.Wo  centered  at  v. 
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FIGURE  5.3 


Figure  5.2  shows  the  results  obtsined  when  the  initial  point,  computed  by  the  con¬ 
tinuation  code  riTCON  (see  [AB83])  with  i/  =  0,  is  a  limit  point  with  respect  to  /i.  Let 
x'  denote  the  dimennonless  radial  deformation  at  the  center.  The  stepsiie  of  the  mapped 
triangles  was  h  =  0.5,  and  the  bounding  hypcrplancs  were  defined  so  as  to  restrict  x‘  to 
the  interval  (1.1, 2.4]  and  i/ to  the  interval  (-0.008,0.006),  respectively.  The  foldline  in 
the  (t/,p)-plane  has  the  shape  pven  in  Figure  5.3.  Figure  5.2,  which  shows  the  manifold 
projected  onto  the  ((',/4)-plaiie,  clearly  shows  a  segment  of  this  foldline  that  includes  the 
local  maximum  and  minimum  points  with  respect  to  i/  at  i/  s  0  and  about  v  =  0.16,  respec¬ 
tively.  The  two-dimensional  simplicial  approximation  algorithm  [/288]  also  captured  these 
two  points,  but  due  to  the  local  nature  of  that  code,  two  runs  were  needed  to  triangulate 
the  manifold  separately  in  the  neighborhood  of  each  of  these  two  points. 


The  third  example  has  beeu  used  in  [A/A/93]  to  test  the  robustness  of  their  two- 
dimensional  code.  The  manifold  is  defined  as  the  subset 

I  (i,v,z)  €  s’  1  F(i,v,r)  =  -z’)  -  i’)-(-e  =  0  }. 

When  <  =  0,  the  cross-section  in  the  (z,r)-plaoe  for  any  fixed  value  of  y  is  the  z-axis  on 
top  of  a  figure-eight,  which  is  symmetric  with  respect  to  the  r-axis  and  intersects  the  z-axis 
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at  2  =  —6,  0,  and  6.  As  <  becomes  positive,  the  solution  set  in  the  (i,i)-plane  breaks 
into  two  components:  one  below  the  z-axis  and  diffeomorphic  to  a  circle,  and  the  other 
the  cross-section  of  an  inverted  trough  above  the  x-axis  and  diffeomorphic  to  the  real  line. 
The  size  of  (  determines  the  width  of  the  trough  opening,  with  larger  values  of  (  giving 
larger  widths.  When  the  variable  y  is  added  to  the  problem  the  result  is  a  two- component 
manifold,  invariant  in  y,  whose  second  component  is  an  inverted  trough  which  is  above 
the  plane  z  =  0  and  whose  opening  is  parallel  to  the  y-axis.  A  narrow  trough  opening  is 
challenging  because  the  solver  may  jump  across  it  without  exploring  the  trough  itself. 


FIGURE  S.4 


Figure  5.4  shows  the  results  of  the  triangulation  algorithm  when  a  =  4.0,  i  =  0.25, 
c  =  5.0  X  10^’,  the  stepsize  is  h  =  0.03,  and  the  initial  point  (z,y,z)  =  (0,0,0.25).  The 
bounding  hyperplanes  were  defined  to  restrict  z  to  the  interval  [-0.15,0.15]  and  y  to  the 
interval  (—0.14,014].  PITMAN  followed  the  curvature  without  falsely  jumping  across  the 
opening  of  the  troiigh. 
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